module compute_single
!DIR$ NOOPTIMIZE
	use compute_gen
	use drawSingle
	implicit none

	real, allocatable	:: rpaths(:,:,:)

	private
	public	:: run_single
	
	contains
	
	function tstat result(stat)
		logical, save	:: init=.true.
		integer, save	:: rind(1000)
		real, save		:: dir(200,10)
		real, allocatable	:: stat(:),v(:)
		integer			:: i,c
		
		if(init) then
			rind=1234872*ranu_v(1000)
			do i=1,size(dir,2)
				call rnnoa(dir(:,i))
			enddo
			init=.false.
		endif
		stat=[log(om2),log(sigma2_F),real(cindF),muU,mt_f]
		c=1
		do i=1,13
			v=[real::]
			select case(i)
			case(1)
!				v=[Xs(:,ir0)]
			case(2)
				v=[f]
			case(3)
				v=[log(kappa2(ir0))]
			case(4)
!				v=[CFs(:,1)]
			case(5)
				v=[rhos(ir0)]
			case(6)
				v=[log(clubs2(1))]
			case(7)
!				v=[CoCFs(:,1)]
			case(8)
				v=[CoCrhos(1)]
			case(9)
				v=[log(CoCs2(1))]
			case(10)
				v=[log(clubs2(1))]
			case(11)
				v=[cindu(ir0)]
			case(12)
				v=[cindCF(1)]
			case(13)
				v=[cindCoC(1)]
			end select
			stat=[stat,v]
		enddo
!		stat=[stat,matmul(transpose(dir(1:size(stat),:)),stat)]
	end function


	
	subroutine draw_once
		call draw_X
		call draw_CF
		call draw_CoCF
		call draw_rho
		call draw_CoCrho
		call draw_kappa2grid
		call draw_clubs2grid
		call draw_CoCs2grid
		call draw_cindu
		call draw_cindCF
		call draw_cindCoC
		if(poosflag) call draw_muU
		call draw_om2
		if(.not. poosflag) then
			call draw_F
			call draw_mt_F
			call draw_fm
			call draw_sigma2_F
			call draw_cindf
		endif
	end subroutine	

	function getpaths result(val)
		real	:: val(npath)
		if(poosflag) then
			val(1:T)=matmul(G(:,0:q0),Us(0:q0,ir0))
			val(T+1:)=fcstUs(1:nh,ir0)+val(T)
		else
			val(1:T)=matmul(G(:,0:q0),f(0:q0)+Us(0:q0,ir0))
			val(T+1:)=fcstf(1:nh)+fcstUs(1:nh,ir0)+val(T)
		endif
	end function

	subroutine iterate
		integer	:: l,leff
		real, allocatable	:: tracestats(:,:)
		
		do l=1,nburnin
			call draw_once
		enddo
!		call mdisp(matmul(G,f))
	call mdisp(matmul(G,f).cvr.matmul(G,Us(:,ir0)))
!	call mdisp(mt_F)
!	call mdisp([sigma2_F,om2,muU])
!stop
!		allocate(tracestats(size(tstat()),ns))
		
		do l=1,nskip*ns
			call draw_once
!			call mdisp(Us(:,5).cvr.Xs(:,5))
!			call mdisp([real(leff),sigma2_F,om2])
			if(mod(l,nskip)==0) then 
				call draw_fcsts
				leff=l/nskip
				rpaths(:,ir0,leff)=getpaths()
!				tracestats(:,leff)=tstat()
				if(mod(leff,500)==0) then
					call mdisp([real(leff),sigma2_F,om2])
!					call mdisp(tracestats(:,leff))
!					call storeML(tracestats(:,1:leff),"ts")
				endif
!				call savestats(leff)
			endif
        enddo
        call printtime
	end subroutine
	
	subroutine run_single
		if (allocated(rpaths)) deallocate(rpaths)
		allocate(rpaths(npath,n,ns))
		do ir0=1,n
!do ir0=1,10
			print *,"now entering single country posterior", ir0
			call mydispnames(names(ir0:ir0))
			call inititerate
			clubid=1
			cocid=1
			if(poosflag) then
				f=postmeanF
				Us(:,ir0)=Xs(:,ir0)-f
			endif
			call iterate
!			call storeML(rpaths(:,ir0,:),"path")
		enddo
		call execinML("clear all")
		if(poosflag) then
			call storeML(rpaths,"Upaths")
			call execinML("save('"//trim(dir)//"single_Upaths_"//trim(cname)//".mat')")
		else
			call storeML(rpaths,"paths")
			call execinML("save('"//trim(dir)//"single_paths_"//trim(cname)//".mat')")
		endif
	end subroutine
	
end module